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We investigate the liquid-vapor interface of the restricted primitive model 
(RPM) for an ionic fluid using a density-functional approximation based on 



> 

in 

m 
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^) , correlation functions of the homogeneous fluid as obtained from the mean- 



spherical approximation (MSA). In the limit of a homogeneous fluid our ap- 



proach yields the well-known MSA (energy) equation of state. The ionic inter- 
facial density profiles, which for the RPM are identical for both species, have 
C , a shape similar to those of simple atomic fluids in that the decay towards the 

bulk values is more rapid on the vapor side than on the liquid side. This is 
the opposite asymmetry of the decay to that found in earlier calculations for 
the RPM based on a square-gradient theory. The width of the interface is, for 
a wide range of temperatures, approximately four times the second moment 
correlation length of the liquid phase. We discuss the magnitude and tem- 
perature dependence of the surface tension, and argue that for temperatures 
near the triple point the ratio of the dimensionless surface tension and critical 
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temperature is much smaller for the RPM than for simple atomic fluids. 
PACS numbers: 68.10.-m, 61.20.Qg 
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I. INTRODUCTION 



In this paper we develop a density-functional theory (DFT) for the properties of the 
liquid- vapor interface of the simplest model of an ionic fluid, namely the restricted primitive 
model (RPM) in which the ions are modelled by equisized hard spheres of equal and opposite 
charge. The RPM serves as a simple model for molten salts and electrolyte solutions. Indeed 
the measured partial structure factors Sij(k), with i,j £ {+,—}, and the thermodynamic 
properties of several molten alkali halides near their melting points can be well described 
by the corresponding quantities of the RPM [Q]. In recent years there has been a revival 
of interest in the properties of the RPM, stemming from efforts to understand the nature 
of criticality in ionic fluids @-§|] • For ionic systems one might suppose that the long-ranged 
Coulomb forces could give rise to critical exponents different from the Ising ones which 
are measured and calculated for atomic and molecular fluids. Some recent experiments on 
certain electrolytes revealed mean-fieldlike behavior or, in some cases, Ising critical regions 
which are several orders of magnitude smaller than in atomic fluids [|J. Since the three- 
dimensional RPM is known to exhibit phase separation into a dense, conducting ionic liquid 
and a very dilute vapor phase, which is also conducting, it is a natural choice for theoretical 
and simulation studies of phase coexistence and criticality. Although attempts to explain the 
experimental observations regarding criticality have so far been unconvincing — estimates 
of the Ginzburg temperatures for the RPM are similar to those for simple (atomic) fluids 
|5||| and the latest Monte Carlo finite size scaling study gives results compatible with 
Ising behavior — this does not mean that the RPM does not warrant further attention. On 
the contrary, because it incorporates the key features of hard core repulsion and Coulomb 
forces it remains the canonical, albeit over-idealized, model for an ionic fluid. 

Here we shift attention away from the bulk and focus on the inhomogeneous situation 
which arises at the planar interface between the coexisting liquid and vapor of the RPM. 
We are interested in the ionic density profiles and the surface tension of such an interface, 
which can be viewed as a crude model for the corresponding liquid-vapor interface of a 
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molten alkali-halide. Compared with the well studied case of simple fluids, modelled by a 
Lennard- Jones potential ||, very little is known about the interface in the RPM. We are 
not aware of any corresponding simulation studies, although there is one early molecular 
dynamics simulation || using the more realistic Born-Mayer-Huggins potential model for 
KC1. 

Theoretical work was pioneered by Telo da Gama et al. JTcJ] who used a gradient expansion 



developed in Ref. |1 !j to investigate the RPM interface. The special symmetry of the RPM 
implies that the density profile of the cations should be the same as that of the anions, 
i.e., there should be local electroneutrality p+(z) = p-(z) throughout the interface. For 
such a "symmetric" situation the gradient expansion of the free energy functional involves 
only gradients of the total density profile p(z) = p+(z) + p~(z); the charge density profile 
q(z) = p+{z) — p-(z) vanishes identically. The coefficients in this expansion involve moments 
of the density-density bulk direct correlation function c£\r; p). If the expansion is truncated 
at the square gradient term, as is usual, the corresponding coefficient must be positive 
if the theory is to yield physical solutions. As was explained in Ref. |TD[ , it is necessary 
to utilize a rather sophisticated, self-consistent theory of the bulk correlation functions, 
the generalized mean-spherical approximation (GMSA) JT^, in order to obtain a positive 
coefficient. The ordinary MSA |T3|JT^1 , which is often successfully employed fl| in studies 



of the bulk RPM, is insufficient as it yields merely the hard sphere (Percus-Yevick) result 
for Cp\r;p) which produces a negative coefficient of the square gradient term. That one 
must supercede the MSA to obtain physical interfaces is, at first sight, quite surprising 
- given its success in bulk. One suspects that this is an artefact of the square-gradient 
approximation. Here we re-investigate the liquid-vapor interface using an alternative DFT, 
which does not utilize the gradient expansion and thereby avoids the need to employ the 
GMSA and the ensuing problems involved with extrapolation into the two-phase region [JXO 



Our approach, which is motivated by the MSA treatment of the bulk free energy, involves 
a local density approximation for the hard-sphere part of the functional and a non-local 
treatment of the remaining (Coulombic) contributions, which is obtained by approximating 



the inhomogeneous pair correlation functions by their homogeneous counterparts. It differs 
from other DFT approaches for ionic fluids [p1j|-|T9|j which have proved successful in the 
primitive model description of electrical double layers at charged hard walls. These theories 
treat the non-hard-sphere part of the functional by means of a second order density expansion 
about the density of a reference fluid, usually taken to be the homogeneous bulk fluid far 
from the substrate. Although this is adequate for many purposes it is problematical when it 
comes to liquid-vapor interfaces or to the adsorption of thick (wetting) films, where two bulk 
phases are involved. Indeed, for the case of an atomic fluid the corresponding second order 
expansion about a homogeneous reference density is known to fail to account for liquid- vapor 



coexistence and is inadequate for wetting problems [20]. We are not aware of attempts to use 



the approaches in Refs. [p^j -19] for the liquid- vapor interface. Nor are we aware of attempts 
to use integral-equation theories for that purpose which have been rather popular in studies 



of the electrical double layer [21,19] but which may be beset by similar problems. 

The present approach does not suffer from these difficulties, i.e., the uniform limit of the 
free energy reduces to that of the MSA for any uniform density pb- On the negative side this 
means that our theory for the interface is prone to the same deficiencies as the bulk MSA, 
namely the failure to incorporate properly the effects of ion pairing, which are especially 
pronounced in the vapor, and the consequent poor estimate of the location of the critical 
point. 

As well as providing a description of the liquid- vapor interface of a near-symmetric alkali 
halide, i.e., one where the ions have nearly the same diameter, the present theory may form 



the basis for a description of wetting phenomena in ionic fluids |22] • The situations one would 
like to consider are: (i) the wetting of a substrate-alkali halide vapor interface by the molten 
salt, (ii) the wetting of the interface between a substrate and, say, the phase dilute in salt 
by the other, salt-rich, phase for an electrolyte that exhibits liquid-liquid phase separation 
such as those considered in Refs. [|J and |23| . Since it is well known that wetting properties 
depend sensitively on the range of the interaction potentials, one might ask if the long-range 
Coulomb interactions lead to qualitatively new features or if, due to screening, ionic fluids 



behave as one-component systems with short range interactions, e.g., Yukawa fluids. 

This paper is arranged as follows. In Sec. [TT] we introduce the approximate DFT for our 
particular problem. Section |TTT| examines the free energy and the second moment correlation 
length which emerge for the bulk fluid, comparing and contrasting the results with those 
of other theories and with exact results for the limit p — > 0. In Sec. [TV] we describe the 
results for the density profile and surface tension of the liquid-vapor interface as a function 
of temperature. These are compared with those of the earlier square-gradient theory and 



with corresponding results for simple atomic fluids. We conclude in Sec. [V| with a summary 
and discussion of our results. 

II. DENSITY-FUNCTIONAL THEORY FOR IONIC FLUIDS 

We study the simplest model for an ionic fluid, the so-called restricted primitive model 
(RPM), which consists of charged hard spheres with equal diameters a for both species. It 
can be considered as a model of a molten salt, but also of an electrolyte solution with the 
solvent treated as a dielectric continuum. The interaction potential is given by 

Wij(r) = w HS (r) + tog(r) (1) 

where 

\ oo, r < a 

w HS (r) = I (2) 
0, r > a 

and 



<(0 = i7 ( r - Q ) ( 3 ) 



where r is the interparticle distance, 9 is the Heaviside step function, i,j e {+, — }, e + = 
— e_ = e is the charge of the particles, and e is the dielectric constant of the solvent. 

The grand-canonical density functional of an inhomogeneous fluid with number densities 
Pi(r) can be written as 



^[{pi(r)}]= F H s[{pt{r)}] + ^Yl j d^d^r' Pi {r) Pj {v) da w^r^gijir, r', {pi(r)}, a) 
-W ^W(r). (4) 

i J 

Here Tus is the free energy functional of the corresponding hard sphere reference fluid, 
r 12 = r — r', Pi is the chemical potential of chemical species i, and denotes the pair 
distribution function of an inhomogeneous system with the density profiles Pi{r) and the 
interaction potential w H (r) + awfAr). The integration over a G [0,1] corresponds to a 



,HS 



straight line in potential space leading from the hard sphere reference potential w to the 
full potential Wij. For the present model this path corresponds to a charging process of the 
ionic liquid. Although Eq. (|4]) is formally exact the pair distribution function is not 
known for the inhomogeneous system. We approximate #y by the corresponding function 
gij(r — r',p,a) of a homogeneous (bulk) liquid evaluated at an appropriate density p, for 
which we choose the mean value of the total densities at the points r and r': 

p = p(r,r')= 1 -J2(Pi(r) + Pi(r'))- (5) 

i 

A similar averaging was used in Ref. [^] in a DFT for a Lennard- Jones fluid in contact 
with a hard wall. However, in that work the local densities were additionally averaged over 
spherical regions around r and r', which is not necessary in the present case because there 
are no pronounced density oscillations at the liquid-vapor interface. Note that the bulk pair 
distribution function g^ depends only on the total density p = p + + p_, which follows from 
the requirement of electroneutrality, p + = p_, in the bulk. (Strictly speaking g i j(r—r f , p, a) is 
not defined for densities which lie within the two-phase coexistence region of the bulk phase 
diagram, but this problem does not arise for the approximate pair distribution functions we 
shall actually use in our calculations below.) 

By writing g^ = 1 + hy the pure Coulombic contribution to the functional 

F C iMr)}] = lT,J rf 3 rrf 3 r' A (r) Pj (r')<(r 12 ) (6) 
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can be separated off. In the bulk, where p + = p_, T° vanishes. Similarly, due to the 
symmetry of the RPM under charge inversion, the density distributions which mini- 
mize Eq. (^) at the free liquid- vapor interface should exhibit local charge neutrality, i.e., 
p + (r) = p_(r) = \p{j). Moreover, the symmetry of the RPM implies that the total pair 

correlation functions hij must satisfy h ++ = h and h + _ = h h , so that Eq. (|]) reduces 

to a functional of the total density p(r): 

Q[{p(r)}]=T HS [{p(r)}}+ 1 - [ d 3 rd 3 r'p(r)p(r') [da—h D (r 12 ,p(v,r'),a) 

1 J Jo er 12 

-p J d 3 rp(r) (7) 

with the difference function h D = (h ++ — h + J)/2 and p = (p + + P-)/2. The hard sphere 
contribution is treated in a local density approximation: 

FhsM*)}] = \\ d 3 r{p(r)[\n(p(r)X 3 ) - 1] + (3f cs (p(r))}. (8) 

A is the thermal de Broglie wavelength and fcs is the non-ideal gas part of the free energy 
density of the bulk hard sphere fluid given by the accurate Carnahan-Starling approximation 

p Arj - 3r/ 2 

fcs(p) = ^ Tr ^r (9) 

with the packing fraction rj = |pa 3 . Clearly the approximation given by Eq. @ does not 
take into account the short-ranged correlations associated with the packing constraints which 
give rise to the oscillatory density profiles encountered for fluids at walls. Nor will it account 
for the weak oscillations that are predicted for low temperatures on the liquid side of the 



liquid-vapor density profile in both simple |2yJ and ionic |28[ fluids. A non-local treatment 
of Ths would be required to describe such oscillatory behavior. 

We employ the bulk correlation function hp given by the analytically solvable mean 



spherical approximation (MSA) [T^JT^ , which provides a reasonable description of the bulk 



structure and thermodynamics of the RPM. Henderson and Smith |2S| have derived an 
explicit expression for ha: 
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Mr, P, T) = J ^T^Y £ ^^ ^l^r ^-^^ - na)r) " J - l((r " " a)r)] 

^ ' n=l ^ ' 

(10) 

where [r/a] is the largest integer smaller than r/a, j n denotes the spherical Bessel function 



of order n and the inverse length T is given by T = (y/l + 2na — l)/(2a) with the inverse 
Debye screening length k = (47r/3e 2 p/e) 1 ^ 2 . The above expression is inconvenient for the 
numerical calculation of hp at large r because in this case many terms have to be evaluated 
and rather large errors may occur due to partial cancellations of the terms. On the other 
hand, a very good approximation for ho for large r can be obtained from the pole analysis of 
the Fourier transform hn(k) as shown by Leote de Carvalho and Evans [p8fl . It provides the 



asymptotes of the form /io(r — > oo,p, T) = A<yp / T ^ exp(— ao(tx)r) for k < k c = 1.228/a and 
hf)(r,p,T) = B<yp ^ exp(— ao(n)r) cos(ai(/«)r + 9(k)) for k > k c with known functions A, B, 
oto, ai, and 6 1 . Using these methods we have derived similar, but slightly more complicated 
expressions for ^-hjj and -§^h D , which are needed to compute the phase diagram and the 
correlation length (see Sec. PH) . These asymptotic approximations have been used in the 
numerical calculations for r/a > 12. 

III. BULK PHASE DIAGRAM AND CORRELATION LENGTH 

For a constant density p(r) = p the density functional yields a Helmholtz free energy 
density 

^ = ^(ln(pA 3 ) - 1) + f C s(p) + \p 2 w Q {p) (11) 

with 

w o(p) — — J cPrr^ 1 J daho{r, p,a). (12) 

Since ho depends only on the dimensionless quantities — and r/a, the charging integration 
with e 2 (a) = ae 2 can be replaced by an integration over (3: 



e 2 r rP 



MP) = J e I d'rr' 1 I dph D {r,p,p). (13) 



This result demonstrates that the free energy given by Eqs. ( |iT|) and (|12l) is equal to that 
obtained by the so-called energy route starting from the MSA pair correlation function. Thus 
the bulk phase diagram is identical to the usual MSA phase diagram as discussed, e.g., in 
Refs. |30|j3l|JT0|j28| j5[1 . The explicit expression for u>o(p) is @ 



\p 2 w {p) = -—i-^^a) 2 + 2na + | - |(1 + 2Kaf 2 ). (14) 



For the reasons discussed in Sec. [TV] below, we also consider the functional which arises when 
we refrain from the integration over a and set a = 1 instead, which amounts to replacing 
the excess free energy T ex = T ' — Ths due to the Coulomb interactions by the corresponding 
internal energy U ex = -^{j3T ex ) (see Eq. (|13|)). We call this approximation scheme MSAl. 
In this case Eq. ( |Tjp is replaced by 

-p W {p) = -— 

2 Anpa 



-p 2 w (p) = —- — 3—^(^0 + (kcl) 2 — nay/1 + 2ko). (15) 



Note that, in contrast to Eq. (|14D , this expression does not reduce to the exact Debye-Hiickel 
result ^p 2 w (p) = — k 3 /(127T/9) in the low-density (k — > 0) limit. However, it differs only by 
a factor of |. 

The liquid-vapor phase coexistence curves, which follow from the free energy by the usual 
double tangent construction, are shown in Fig. [I] for the approximations MSA and MSAl, 



together with Gibbs ensemble Monte Carlo simulation results [p2||33| . The results are given 
in terms of the reduced density p* = pa 3 and the reduced temperature T* = k B Tea/e 2 . The 
phase diagram of the RPM is characterized by small values of the reduced critical density 
and temperature (for comparison the Lennard- Jones fluid has T* ~ 1.3 and p* ~ 0.3 in 
Lennard- Jones reduced units) and by a strong asymmetry of the coexistence curve. At low 
temperatures the vapor phase becomes extremely dilute, e.g., at T* = 0.04 ~ T*/2 one has 
p* ~ 10~ 7 within the MSA. It is well known that the MSA overestimates the critical tem- 
perature and grossly underestimates the critical density. The cruder approximation MSAl 
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predicts an even higher critical temperature (MSA: T* = 0.0786, MSAl: T* = 0.0846) 
and smaller critical density (MSA: p* = 0.0145, MSAl: p* = 0.0086). Simulation results 
for the coexistence curve have undergone substantial revision during recent years (see, e.g., 
Refs. [@,§J); the latest estimate of the critical point by Caillol et al. J7J, using a mixed- field 
finite size analysis of Monte Carlo data, gives T* = 0.0488 ± 0.0002 and p* = 0.080 ± 0.005. 
Clearly quantitatively reliable results cannot be expected from either MSA or MSAl. Both 
fail to take proper account of the effects of ion pairing which is known to be very strong 
in the RPM at low densities and which is believed to influence strongly the location of the 
critical point pQ|p|fl. Since the differences between the coexistence curves obtained from 
MSA and MSAl are relatively minor we do not expect to introduce significant additional 
errors by using MSAl. We note that the main contributions to the integration over a in 
Eq. (|j"2"D stem from the region near a — 1. 

The chemical potential p at coexistence is, for a given temperature, 

d F(p) 



d F(p) 

^0 — TT 



(16) 

Pv 



dp V 

where pi and p v are the densities of liquid and vapor, respectively. Knowledge of po(T) is a 
prerequisite for the subsequent analysis of the liquid- vapor interface. 

We now turn attention to the bulk correlation length that follows from the present theory. 
As we shall see in the next section, the width of the interface is governed by the correlation 
length of the bulk liquid. Every density functional ansatz defines a direct correlation func- 
tion via twofold functional differentiation of the intrinsic Helmholtz free energy functional 
jF[{p,j(r)}] with respect to density pi 



^( r y, Wr ) } ) = H9 f2^ + Ml^. (17) 

6pi{T)6p j {T f ) pi(r) 

From the bulk limit of this function a pair distribution function hij can then be obtained via 
the Ornstein-Zernike equation. The results will be different from the used for the con- 
struction of the functional, in our case those of the MSA. This is an important advantage 
of the density-functional approach. We require the correlation function for fluctuations of 
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the total density hf^ = (h ++ + h+-)/2 to exhibit the conventional Ornstein- Zernike behav- 
ior, i.e., hf^ should become long-ranged near the critical point, with an exponential decay 
described by a diverging correlation length. Such a behavior reflects the net inter-particle 
interaction which is necessary to obtain liquid-vapor coexistence. However, as remarked in 
the Introduction, the original MSA result for h p is simply that of hard-spheres (in the PY 
approximation), so there is no manifestation of attraction in this particular combination of 
the correlation functions. We define the second-moment correlation length £ for fluctuations 
of the total density by the expansion of the structure factor: 

S(k, p) = l + phf\k, p) = S(0, p)/(l + i 2 k 2 + •••), (18) 

where k is the wave number. This correlation length can be obtained rather easily if one 
again assumes p+(r) = p_(r) = |p(r) and considers T = Q + p J d 3 rp(r) (see Eq. (|^)) as a 
functional of the total density p(r). Using the bulk limit c p of the function 



£ is given by 



2 = 1 fd 3 rr 2 c p b) (r,p) 
6l/p- Jd 3 rc p b \r, P y 



which follows straightforwardly from the definition given by Eq. pq ) and the Ornstein- 
Zernike equation h p b \k,p) = c p (k,p)/(l — pc p b \k,p)). Within the MSAl approximation 
one finds 

/rtf (r,p) = -5(v)^f cs (p) - ^-[h D (r, p) + p^-h D (r, p) + \p 2 ^h D (r, p)} 

-8(r) J dh'^[p-^h D (r',p) + \p^h D (r',p)\. (21) 



Note that the equivalent result for a simple atomic fluid was derived in Ref. ||34j| . The re- 
sulting correlation lengths of the coexisting liquid and vapor phase are plotted in Fig. 0. 
In accordance with the mean-field character of the present approach they diverge near the 
critical point proportional to (T c — T) 1//2 . £ also diverges upon approaching the spinodals, 
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which are determined by 9 gjf^ = 0, or, equivalently, by the vanishing of the denominator 
in Eq. (^0|). Upon lowering the temperature the correlation length in the liquid phase de- 
creases to about one particle diameter, whereas that of the vapor phase increases again as 
a consequence of the reduced screening at very low densities. 

It is instructive to examine the limiting behavior as p — > 0. The exact low density 



behavior of the second moment correlation length is believed to be |35 

1/4 



H(S i i+ °(" 1/2 )i- < 22 > 



a result which also follows from the HNC approximation p6|j35| . If one assumes that for 
low densities the MSA correlation function ho tends to — ^re~ Kr all the algebra can be 
performed analytically and the MSAl yields in leading order 

which is in accordance with our numerical results. Thus, for p — > the MSAl correlation 
length diverges with the correct power law but the amplitude is too large by a factor of \fl. 
We can also carry through the algebra using the original functional Eq. ([/]) (including the 
a integration). In this case we obtain 

^ 4 \36nepJ V ; 

for p — > 0. Once again the power law is correct but the amplitude is again too large. 
These results suggest that the functionals capture the essential features of density-density 
correlations in the p — > limit. This is a non-trivial observation. Recall that GMSA, which 
assumes a particular (single Yukawa) form for the non-Coulomb part of the direct correlation 
functions Cjj(r) and then enforces consistency among three routes to bulk thermodynamic 



functions, yields a finite value as p — > p5[. Thus our present MSAl and MSA results, 



Eqs. ( p3]) and (p4]), improve upon the GMSA in the low-density limit. Finally we should 
emphasize that the results described here refer to the second moment correlation length £ 
as defined by Eqs. (pl|) or (p0|). The so-called true correlation length ^ that determines the 
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ultimate exponential decay of the density- density correlation function h^\r, p) is determined 
by the poles of the Fourier transform hf\k, p) ||28|| . In general ^ differs from £; they become 
identical (within mean-field theory) in the vicinity of the critical point. Calculating the poles 
from the present DFT approach is not straightforward as the Fourier transform of Eq. ( |2~T|) 
cannot be performed analytically, so we cannot easily compute . However, we are able 
to compare our results for £ with Coo,gmsa obtained in Ref. P8[ ]. For subcritical liquid-like 
states which turn out to be relevant for the interface, the correlation lengths £msa an d £,msai 
both differ by less than a factor of two from ^oo,gmsa- In the very low density vapor phase 
the exact limiting behavior is ^ ~ 1/(2k) |36|j35|] which diverges as (T/p) 1 / 2 for p — ► 0, i.e., 
faster than the second moment correlation length given by Eq. (0). 

In summary we conclude that the present density-functional theory provides a reasonable 
description of the long-wavelength (small k) behavior of density-density correlations in the 
bulk RPM. This is a necessary prerequisite for a reliable treatment of the liquid- vapor 
interface. 



IV. LIQUID- VAPOR INTERFACE 
A. Density profiles 

The average density profile of the liquid-vapor interface varies only in the direction 
normal to the interface which we take as the z direction. In the following we assume that 
the bulk liquid density pi is obtained for z — > — 00 and the bulk vapor density p v for z — > 00. 
The application of the density-functional theory outlined in Sec. || to this situation allows 
one to carry out the lateral integrations in Eq. ([/]) yielding the following grand-canonical 
functional per surface area A for the total density p(z) = 2p_(z) = 2p + (z): 

^tt[{p(z)}] = ^Fhs[{p(z)}\ + \ J dzdz' p{z)p{z')w{z - z',p(z,z')) - p J dz p{z) (25) 

with 
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w(zi2,p) = / dr dahi>(r, p,a). (26) 

6 J\Z12\ Jo 



.,2 poo pi 

\ dr 

'\Z12\ JO 

For densities not too large the function hp exhibits its slowest decay as function of r 
(lowest value of a (n), the imaginary part of the leading pole) for small values of the charge, 
i.e., for small a (see Fig. 4(b) in Ref. ||28|| ). Therefore, in contrast to the bulk free energy, 
the large distance behavior of w(zi2,p) is dominated by the contributions from small a. In 
this region the correlation function is expected to take on the Debye-Hiickel form 

Mr, p, a) ~ -Pf-ae-*^ (27) 
er 

which also follows from the small coupling expansion of the decay length and of the expression 
for the amplitude obtained from the pole analysis []28| mentioned in Sec. fll. From Eqs. ( |26|) 
and (1271) one finds 



^-w(z — oo, p) ~ [ da ae~ K ^ z (28) 

UZ 6 Z Jq 



and hence 



. . 67r/?e 4 

w (z^oo :P )^-^-^. (29) 

Thus w, which corresponds to the effective interaction between two fluid layers a distance 
z apart, is predicted to decay algebraically with the same exponent as for a Lennard- Jones 
fluid, where the interatomic potential decays as r~ 6 . This unexpected result implies an 
algebraic decay for the density profiles, too. More specifically one would expect dp/dz ~ 
— \z\~ A as z — > ±oo. This behavior can be traced back to the fact that within the MSA 
density-functional theory the bulk direct correlation function cf\r;p), which is given by 
Eq. ( pl~D with hz>(r, p) replaced by L dahn(r } p,a), also decays algebraically. Using an 
argument similar to the one above one finds cf (r,p) ~ r~ 6 as r —>■ oo, which is the same 
as for a Lennard- Jones fluid. This would, in turn, imply that hp(r;p) as obtained from 
Cp (r, p) via the Ornstein-Zernike equation should also decay as r~ 6 . In other words the 
theory generates effective potentials for density-density correlations which mimic those of 
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genuine dispersion forces. Whilst this observation is certainly intriguing it is likely to be 
an artefact of the present MSA scheme. The bulk pair correlation function hp (r\ p) should 
decay exponentially for the RPM [^[^J. Moreover, general arguments for the decay of 
density profiles near walls P7| , |28| , |3"7| show that the ultimate decay of p(z) should mimic that 



of hp. For example p\ — p(z) should decay as exp(z/^oo) as z — > — oo, where £00 is the true 
correlation length of the bulk liquid. If, on the other hand, we ignore the a integration and 



set a = 1 (see also the discussion in Sec. |TJ) then — dw ^ z ' p ^ ~ h D (z, p) and there is no longer 



an algebraic decay of the direct correlation function and of the profile. Henceforward we 
describe results based on both MSA and MSAl. 

By functional differentiation of Eq. ( p5|) one derives the Euler-Lagrange equation 

logp(z)A 3 + f cs (p(z)) = po-p(z) (30) 

with 



p(z) = / dz p{z 



w(z - z', p(z, z)) + -p(z)w'(z - z', p(z, z')) 



(31) 



(The primes on fcs and w denote derivatives with respect to the density.) In contrast to 
similar theories for Lennard- Jones or dipolar fluids the kernel in Eq. ([IT]) depends 



on z and z 1 separately due to the density dependence of the approximated pair distribution 
function. However, in the present case one can still show that if a profile p{z) solves Eq. ( |30|) 
any shifted profile p{z — z ) will be a solution as well, because the position of the interface 
is not fixed by the boundary conditions. In order to select a unique profile we have imposed 
the additional requirement p{z = 0) = (pi — p v )/2 and have shifted the profiles accordingly 
during the numerical solution. 

Equation ( |30"D was solved using the usual Picard iteration scheme. In each iteration (n) 
the function p(z) is calculated from the density profile P( n )(z) in the region z G [— L/2, L/2] 
and a new p new {z) is obtained by numerical inversion of Eq. ([30]) for each value of z. In 
order to achieve convergence the profile for the next iteration is calculated according to the 
mixing rule 
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p (n+ l) = LOp new + (1 - U})p(n) (32) 

where typically u = 0.2. The process is repeated until max 2 \p( n+ \)(z) — p( n )(z)\ is smaller 
than a prescribed accuracy. Since the functions w(zi 2 ,p) and w'(zi2,p) have discontinuous 
derivatives at zu = ±a the integration in Eq. (|31~|) is divided into the three subintervals 
[—L/2, 2 — a], [z — a, z + a], and [z + a, L/2]. The Milne rule, which effectively interpolates 
the integrand piecewise by third order polynomials, is used for all integrations. Asymptotic 
contributions to p(z) from the regions \z'\ > L/2 are calculated by replacing p(z') by the 
bulk limits p{z! < —L/2) = pi and p(z' > L/2) = p v and integrating numerically with 
w(\z - z'\ > 20a) = 0. 

The functions w and w' are given by integrals over the functions ho and h' D whose 
evaluation already requires significant numerical efforts. Since w(zu,p) and w'(z\2,p) have 
to be evaluated many times during the iteration scheme the algorithm can be accelerated 
considerably by using two-dimensional spline interpolations for these functions, which require 
their evaluation only at, e.g., 100 x 100 grid points before the main algorithm starts. 

In Fig. [| we plot the density profiles obtained from MSAl for a series of different tem- 
peratures. Upon increasing the temperature towards the critical point (T* = 0.08465) the 
interface broadens and the density difference pi — p v between the phases decreases. These 
obvious features can be incorporated by introducing a scaling function p sc i(z/£,T) with 
p scl (±oo,T) = =Fl: 

p(z, t ) = \{pi + Pv) + ^(Pi ~ Pv)Psd(z/£,T). (33) 

For T —>■ T c the scaling function should reduce to a universal function of the single variable 
z/C,, which — within mean-field theory — is given by 

p scl (z/£, T^T C ) = - tanh ± (34) 

Such a plot is given in Fig. |4] using the second moment correlation length £z(T) for the bulk 
liquid phase, obtained from Eqs. and (pT|). This demonstrates that the scaling behavior 
remains valid even relatively far outside the critical region. In order to make comparison 
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with the behavior of simple atomic fluids we present in Fig. [5] the density profiles of a 
Lennard- Jones fluid scaled in the same manner. These results have been obtained using 
the density-functional ansatz of Frodl and Dietrich . At low temperatures the profiles of 
the ionic fluid deviate weakly from antisymmetry, approaching the bulk limit faster on the 
vapor side than on the liquid side. The same trend can be observed in the Lennard- Jones 
system, where it is even more pronounced: compare profiles corresponding to the same 
reduced temperature t — 1 — T/T c . In the case of a Lennard- Jones fluid the deviation from 
antisymmetry is attributed |4(| to the fact that fa < fa at low temperatures. The situation is 



different for the ionic case where, away from T c , fa ^> fa (see Fig. |2|). This would imply that 
within a square gradient theory there should be a slower decay of the profiles on the vapor 



side. The results of Ref. [ 10 1 are in accordance with this expectation. However, compared 
with the square gradient theory our present non-local theory clearly predicts the opposite 
asymmetry: the decay on the vapor side is faster than on the liquid side, although fa ^> fa. In 
order to understand these features it is important to distinguish between the intermediate 
range decay that is apparent from Figs. and the ultimate asymptotic decay into the 
bulk. As mentioned previously, for the ionic fluid the ultimate decay should be exponential 
with a decay length equal to the true correlation length of the bulk. We attempted to 
analyze the numerical results for the tails of our profiles but this is rather difficult because 
the variations are very small. When curves of the form Cexp(— |z|/£) or (C/z)exp(—\z\/fa) 
are fitted to the tails of the MSAl profiles one finds larger values for £ on the vapor side, as 
expected, but the actual values depend very strongly on the z interval used for the fitting. 
Finally we note that the density profiles in a Lennard- Jones fluid approach, as a function of 
temperature, the universal scaling function monotonically from below on the vapor side and 
from above on the liquid side (see Fig. |5|). On the other hand in the ionic fluid the density 
profile on the vapor side is below the universal profile at low temperature but is above it 
for large z close and T to T c (see Fig. |] for t = 1 — T/T c = 0.0313). We interpret this as a 
consequence of the facts that fa > fa and that the vapor becomes more important near T c 
due to its increasing density. At the same reduced temperature the scaled interfacial profiles 



in a Lennard- Jones fluid are steeper than those of the ionic fluid. 
We define the interfacial width 5 as 



(35) 



2=0 

For the tanh profile in Eq. ( |3~3D one has 5 = 4£, so that in the critical region 5/4 should 
diverge in the same way as the correlation length. Remarkably, Fig. |2| shows that at all 
temperatures this width is determined by the liquid correlation length £j. At T* = 0.03, 
which is a little higher than the estimated triple point temperature ||33|| , S/a ~ 3.0. This 
value is slightly higher than the estimates of the equivalent ratio for simple atomic fluids 



near their triple points ||. Also shown in Fig. are the results of Ref. [10] which exhibit 



a similar variation with T but at low T the widths are about a factor of two smaller. The 



'10-90' interfacial thickness used in Ref. [1C] is defined as the distance over which the total 
density changes from p v + 0.9(pi — p v ) to p v +0.1(pi — p v ). For the tanh profile it gives a width 
that is larger by a factor artanh(0.8) = 1.0986 than the quantity 5 defined in Eq. fl35|). The 



data taken from Ref. [nj have been divided by this factor. For the actual profiles shown in 
Fig. ^| the two definitions for the width yield values which differ by about 10%, too. 

The density profiles obtained from the MSA (including the a integration) are, for a given 
value of t — 1 — T/T c , similar to those obtained from MSA1. When plotted in terms of scaled 
variables, Eq. (^), the profiles from the two theories can hardly be distinguished on the scale 
of Fig. |4j. Moreover, identifying the ultimate algebraic decay from the numerical results for 
the tails of the profiles was not feasible. The interfacial widths 5 calculated from the MSA 
profiles are close to 4£j, with obtained within the MSA approach. These widths are in 
turn close to those from MSAl, provided we make comparison at the same value of t. 

B. Surface tension 

In this subsection we present results of calculations of the liquid-vapor surface tension 7. 
First we derive an expression for 7 in terms of the density profile p(z) which is convenient 
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for numerical work. Since it is necessary to cut off the system at z — ±L/2 this generates 
artificial fluid-vacuum surface tensions 7i,« ac and 7„ )Uac which must be subtracted in order to 
obtain the liquid-vapor surface tension. Therefore 7 is given by 



7 = lim 

L^oo 



7i,fctc 7t> 



(36) 



A A 

where Q(L) is the grand-canonical functional evaluated for a finite system of size L. It is 
convenient to express the bulk free energy (Eq. (|TT|)) as a double integral over z and z' 
analogous to Eq. (p5|). A straightforward calculation yields 



L/2 



dz 



L/2 



dz' w(z — z', p) = — wo(p) — 2 / dzt(z,p) 



L/2 



(37) 



with t(z,p) = dyw(y, p) and, from Eqs. ( |TB| ) and 

/■oo 

w o(p) = 2 / dyw(y,p) = 2t(0,p). 



(38) 



The surface tensions with the vacuum follow from assuming constant density profiles p(z) 
Pb, with b = v,l, which leads to (there are two equal interfaces) 



lb,vac = ■= lim -p b 
Z -L— >oo Z 



L/2 



L/2 



L/2 



L/2 



dz' w(z - z', p b ) - Lw (pb) 



(39) 



Using Eq. ( [371) and the symmetry w(—y) = w(y) gives 



76, 



■■■>c = TiPb r lim 

Z L— >oo 



oPf, r lim 

Z L— >oo 



L/2 



dz' w(z — z' , pb) — 2 / dzt(z,pb) 



-L/2 



L/2 



L/2 



dz(t(z,p b ) -t(L/2 + z,p b )) -2 / dzt(z,p b ) 



L/2 



(40) 



/ dzt(z,p b ). 



Thus the final expression for the surface tension follows by inserting Eqs. (|TTD, (p7|), and 
(M) into Eq. 



7 



J /"OO 



dz[f H s(p(z)) - fHs(pSK(z))] - fio dz[p(z) - Psk(z)} 

J —00 

POO 

dz / [p(z)p(z')w(z - z',p(z, z')) - plw(z- z',p v )] 
Jo 
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(41) 







+ - I dz I dz' [p(z)p(z')w{z - z',p(z,z')) - pfw{z - z',pi)] 



oo J — oo 
/-oo 



roo 2 /*oo 

+ / dz I dz'p(z)p(z')w(z - z',p(z,z')) - -pi / dzt(z,p v ) - -pf / dzt(z,pi) 

J-oo JO 1 JO 1 Jo 

with the sharp-kink profile 

, Ph z < 

Psk(z) = r ( 42 ) 
z > 

and fnsip) — ^ [ln(pA 3 ) — 1] + fcs(p)- All integrals occuring in Eq. (|4l|) are convergent. 

Since the surface tension is the surface excess grand potential per unit area, an alternative 
starting point is the formula 

/oo 
dz(u{z)+p) (43) 
-oo 

with the grand-canonical potential density 

1 f°° 

w(z) = fns(z) - Pop(z) + -p(z) / dz'p(z')w(z - z, p(z, z')) (44) 

J -oo 

and the vapor pressure 

V = -fHs(Pb) + P-oPb - ^p 2 b w (p b ), b = l,v. (45) 



It is straightforward to show that Eq. (|43D also leads to Eq. ([41]) . 

The surface tensions obtained from the MSAl and MSA theories are displayed in Fig. |6| 
as a function of T/T c and are compared with the results of the square- gradient theory ||10|| . 
Since MSAl has a higher critical temperature T c than MSA we plot the ratio 7 = 7a 2 / (fcsT c ) 
rather than the reduced tension 7* = 7ea 3 /e 2 , which is displayed in Ref. |10| . The MSA 



results for 7 are larger than the values predicted by MSAl which in turn are larger than 
those obtained from the square-gradient theory, with the differences becoming larger at low 
temperatures. Near T c all three theories yield the standard mean-field critical behavior, i.e., 
7 ~ (1 — T/T c ) 3//2 . At lower temperatures the square-gradient theory gives a nearly linear 
variation of 7 with T/T c , whereas in MSAl the (1 — T/T c ) 3 / 2 behavior appears to persist 
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further from T c , which is consistent with the earlier observation that the density profiles 
follow the scaling behavior even for temperatures well removed from T c . 

The present MSAl results predict 7(T/T C = 0.6) ~ 0.09. This estimate should be con- 
trasted with the corresponding result for a Lennard- Jones fluid where simulations and the- 
ories yield 7(T/T C = 0.6) ~ 0.6 ||. The physical origin of such a large difference is not 
obvious and it is important to ask whether it is an artefact of the present approximations. 
Although our theories overestimate T c they should make a compensating overestimation of 
the surface tension 7. It is conceivable that since MSAl and MSA grossly underestimate, for 
a given value of T/T c , the simulation results for the liquid densities at coexistence (see Fig. [I]) 
they also underestimate the magnitude of 7. On the other hand, if one takes experimental 
data for the surface tension of molten alkali halides near the melting points and makes some 
reasonable estimates of the (average) diameter a the values of the resulting reduced tension 



7* [10] are similar to those obtained by extrapolating the present results to the appropriate 
reduced temperatures T*. This would suggest that the present theories yield at least the 
correct order of magnitude for the surface tenstion 7. 



V. SUMMARY AND DISCUSSION 

We have developed two, closely related, density-functional theories for the liquid-vapor 
interface of the RPM. Both the MSA and the MSAl are based on approximating the in- 
homogeneous pair correlation function by that of a homogeneous bulk liquid at some mean 
density p. The MSAl involves an additional approximation, i.e., the integration over the 
coupling constant (charge) is ignored. The results which have emerged from this analysis 
can be summarized as follows: 

1. MSA and MSAl yield similar bulk phase diagrams, with the latter overestimating the 
critical temperature T c and underestimating the critical density p c to a further extent 
than MSA (Fig. H). 
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2. The two theories yield similar second moment correlation lengths for the density- 
density correlations in the bulk fluid (Fig. |2|). They exhibit the correct power law depen- 
dences in the low-density limit but overestimate their numerical prefactors (Eqs. (|2"2"|)- 

3. When suitably scaled to take into account the differences in the bulk phase diagram 
the total density profiles are similar in both theories. Both predict that for low temper- 
atures and on intermediate length scales the profile approaches its bulk limits faster 
on the vapor side than on the liquid side (Figs. § and §). This is opposite to what is 
found in the square-gradient theory of Ref . |I| . However, the ultimate decay into bulk 



should be determined by the bulk correlation length which is larger on the vapor side. 

4. The interfacial width 5, defined by Eq. fl3"5|), is close to four times the liquid correlation 
length over the whole temperature range. It is larger but has a similar variation with 
temperature as that of a simple liquid (Fig. 0). 

5. For T — > T c the density profiles approach a universal scaling form. In contrast to simple 
fluids this approach is not monotonic as a function of temperature (Figs. | and |5|). 

6. The surface tensions calculated 7 from the present theories are larger than those ob- 



tained in Ref. |TI| and the temperature dependence is not as linear (Fig. 



7. All theories for the RPM find that for T/T c ~ 0.6 the scaled surface tension 7 = 7*/T* 
is considerably smaller than the corresponding ratio for simple atomic liquids. We 
surmise that this is a genuine feature of ionic fluids. 

Certain aspects of our theory warrant further discussion. The basis of our approach is the 
approximation for the inhomogeneous pair distribution function <7y(r, r', {pj(r)}, a) which 
is the key input for the theory. This has two intriguing consequences for the asymptotics of 
the correlation functions generated from the density functional. The first concerns the low- 
density behavior of the second moment correlation length £, referred to above. That MSA 
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and MSAl do not yield the exact limiting behavior for £ reflects the fact that even at low 
densities our approximation for does not capture the proper long wavelength variation, 
i.e., the correct coefficient of k 2 in the structure factor (Eq. (0)). A possible improvement 
upon the present approximation for ^ (r, r', {/^(r)}, a) may be given by the "mean density 
approximation" (MDA) used in Refs. |33J and in which the inhomogeneous pair distribu- 
tion function is expanded around its homogeneous limit up to second order in the deviation 
p(r) — pb- But, as for DFTs based on a density expansion of the excess free energy, this 
approach cannot be applied straightforwardly to the liquid-vapor interface due to the lack 
of a unique bulk density that could serve as a starting point for the expansion. The second 
intriguing consequence is that the MSA, through its integration over a in function space 
(Eq. (H)), yields algebraically decaying bulk correlation functions and liquid- vapor density 
profiles whereas the correct ultimate decay should be exponential in both cases. Once again 
this failing must be attributed to the simple approximation employed for the inhomogeneous 
pair correlation function. However, it is not clear what modifications to the theory should 
be introduced — other than simply omitting the a integration as was done in MSAl — in 
order to eliminate this feature. 

We should also emphasize that the present DFT of the liquid- vapor interface is a mean- 
field treatment in that it does not incorporate the effects of either bulk critical fluctuations 
or of interfacial capillary- wave-like fluctuations. In keeping with treatments of simple atomic 
fluids we argue that the density profiles and surface tensions we calculate from our theory 
are those of the bare or intrinsic interface. In order to estimate the additional broadening of 
the density profile arising from the interfacial fluctuations one could perform the standard 
Gaussian unfreezing of these fluctuations on the intrinsic interface given by the present 
theory [fjSpBf . The "stiffness" of the interface is determined by the dimensionless ratio to = 
fc^T/(47T7^ 2 ), i.e., the larger the surface tension 7 the smaller is the interfacial broadening. 
Since we predict that in the RPM the ratio 7 = 7a 2 /(fc^T) is much smaller than for an 
atomic fluid at the same reduced temperature T/T c , this implies that the broadening due 
to fluctuations is significantly more pronounced for the RPM. 
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There are two other interesting problems to which the present theory could be applied. 
The first is the interface between an ionic fluid and a charged wall. Although this has been in- 
vestigated using several theoretical techniques [^T],[n| the advantage of the present approach 



is that it incorporates two coexisting phases and thus can describe wetting phenomena. The 
second is the liquid-vapor interface of an ionic fluid in which cations and anions have un- 
equal diameters. These more realistic systems have only been tackled within the context of 



the gradient expansion [43]. In both cases electrical double layers will form giving a local 
violation of electroneutrality. However, the bulk correlation functions used as input to the 
DFT are only available for neutral systems. Nevertheless, one might still hope to capture the 
essential features by employing Eq. @ for the mapping to the bulk system and by taking 
account of the non-zero local charge density only via the Coulomb contribution to the free 
energy, Eq. (^]). A similar assumption has proved successful for the electrolyte- wall interface 
in the density-functional theories of Refs. [fL^-18 . 
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FIGURES 

FIG. 1. The phase diagram of the restricted primitive model in the temperature-density plane as 
given by the approximation schemes MSA and MSAl discussed in the main text. The open symbols 
denote the results of Gibbs ensemble Monte Carlo simulations |3^]. The solid circle denotes the 
most recent simulation estimate of the critical point M. The vapor-liquid-bcc-solid triple point is 



estimated to lie near T t * ~ 0.025 and pi ~ 0.5 33]. The dotted line connecting the Monte-Carlo 
data is a guide to the eye. 

FIG. 2. The second moment correlation lengths and for density correlations in the vapor 
and liquid phases at two-phase coexistence obtained using the MSAl approximation. Also shown 
is the width 5 (defined by Eq. (|35|)) of the liquid- vapor interface obtained from the present den- 



sity-functional theory and from the square gradient theory of Ref. [10]. All lengths diverge at the 
critical temperature T c . The vapor correlation length increases at low temperatures because of the 
reduced screening in the very dilute vapor phase. The inter facial width 5 is very close to 4£; over 
the whole temperature range. At the same reduced temperature the interfacial width predicted by 
the present theory is larger than that obtained from the square-gradient theory which in turn is 
close to the result for a Lennard-Jones fluid as obtained from a DFT similar to the present one 
& 



FIG. 3. The total density profile of the liquid-vapor interface of the RPM at different temper- 
atures, as obtained from MSAl. 

FIG. 4. Density profiles, as obtained from MSAl, plotted in scaled form (see Eq. (j33|)) as a 
function of z/&. The parameter t = 1 — T/T C measures the deviation from the critical temperature. 
For t — > the scaled profiles approach the universal scaling function — tanh(z/(2^)). 

FIG. 5. The same plot as Fig. Qbut for a Lennard-Jones fluid. The results were obtained using 



the density- functional theory of Frodl and Dietrich [39|. Note that the approach to the scaling limit 



is slower and the profiles are more asymmetric than for the ionic fluid. 
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FIG. 6. Comparison of the liquid-vapor surface tension 7 of the RPM as given by the present 



theories MSA and MSA1 and by the square gradient theory of Ref. [1C]. The data have been 
reduced using the critical temperature T c appropriate to each theory. 
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